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Despite the significant progress made in recent years, the computation of the complete set of 
elementary flux modes of large or even genome-scale metabolic networks is still impossible. We 
introduce a novel approach to speed up the calculation of elementary flux modes by including tran- 
scriptional regulatory information into the analysis of metabolic network. Taking into account gene 
regulation dramatically reduces the solution space and allows the presented algorithm to constantly 
eliminate biologically infeasible modes at an early stage of the computation procedure. Thereby, 
the computational costs, such as runtime, memory usage and disk space are considerably reduced. 
Consequently, using the presented mode elimination algorithm pushes the size of metabolic networks 
that can be studied by elementary flux modes to new limits. 



INTRODUCTION 

Elementary flux modes (EFM) are indivisible sets of re- 
actions that represent biologically meaningful pathways 
P3 H] under steady state condition. Removing only a sin- 
gle reaction of an EFM results in the extinction of the 
entire pathway. Consequently, EFMs can be used to de- 
compose metabolic networks mathematically and investi- 
gate them unbiasedly. For that reason EFMs have gained 
increasing attention in the field of metabolic engineering 
in recent years. However, the computation of EFMs is 
of combinatorial complexity [3 . Hence, the computa- 
tional costs for calculating EFMs increase sharply with 
the size of the analyzed network. The calculation of all 
EFMs of small networks (up to 50 reactions) is straight- 
forward and simple. However, despite the major progress 
made recently [4 - 6 the computation of the complete set 
of EFMs of large or even genome-scale networks is still 
impossible. There is a number of tools specifically de- 
signed to calculate the complete set of EFMs as perfor- 
mant as possible, such as Metatool [7], CellNet Analyzer 
[5] and efmtool [5]. To our best knowledge one of the 
fastest program currently available is efmtool by Marco 
Terzer which is written in the multi-platform program- 
ming language Java, supports multi-threading, is pub- 
lished under the open source software license Simplified 
BSD Style License j5] , and can be downloaded from |10j . 

In the presented work we introduce a novel approach 
to speed up the computation of the complete set of bio- 
logically feasible EFMs by taking into account the gene 
regulatory information of the investigated metabolic net- 
work. Transcriptional regulatory networks (TRN) are 
typically provided as a boolean rule set, e.g. [TT]. These 
rules exclude many of the mathematically possible EFMs 
for biological reasons. We implemented our algorithm by 
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extending efmtool, thereby, exploiting the full power and 
advantage of open source software. By utilizing a specific 
feature of the binary approach [3] which was applied in 
efmtool, the elimination of biologically infeasible modes 
can be done constantly and at an early stage of the EFM 
computation process. Thereby, a significant reduction of 
the computational costs, such as execution time, memory 
consumption and harddisk space, is achieved. 

METHODS 
Binary approach 

Modern EFM computation programs, such as efm- 
tool, use a binary approach [4] of the double description 
method |12j . In the following we briefly review this bi- 
nary approach. We will introduce our modifications for 
the inclusion of transcriptional regulation in the next sec- 
tion. 

The binary approach is characterized by splitting each 
mode into a binary part and a numerical part. The bi- 
nary part of a mode contains only a single bit for each 
reaction, where T' means that the reaction carries a flux 
and '0' stands for a reaction not carrying a flux. While it- 
erating through the binary algorithm, the numerical part 
of each mode is successively converted into the binary 
representation. The iteration procedure terminates when 
each mode has been completely transformed to its binary 
form. In a final post-iteration step the computed binary 
modes are converted back to their numerical forms. The 
numerical representation of a mode gives the exact stoi- 
chiometric amount of each involved reaction that partic- 
ipates in the mode. 

We demonstrate the general principals of the binary 
approach by the simple example network shown in Figure 
[T] For the sake of clarity the network will not be com- 
pressed in order to keep all originally specified reactions 
and metabolites of the network. In a 'real-life' computa- 
tion several compression strategies would be applied first 
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FIG. 1. Example network consisting of 12 metabolites (rect- 
angles) and 11 reactions (diamonds). Only one of the 11 re- 
actions is reversible (R7r). 
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TABLE I. Extended stoichiometric matrix S ex t of the exam- 
ple network shown in Figure [T] after splitting the reversible 
reaction R7r into the two irreversible reactions i?7/and R7b. 



in order to combine and remove topographically redun- 
dant reactions and metabolites [3] . The example network 
consists of 11 reactions and 12 metabolites. Only reac- 
tion R7r is reversible. The stoichiometric matrix S of the 
example network is shown in Table 1 of the supplemen- 
tary data section. The external metabolites do not obey 
the steady state condition and, thus, are irrelevant for 
the calculation of the EFMs. 

First, the reversible reaction R7r is split into a for- 
ward and a backward irreversible reaction. This is done 
by negating the column of the reversible reaction and ap- 
pending the newly created column right after the original 
one. Tabic U shows the extended stoichiometric matrix 
S ex t that only contains irreversible reactions. 

The main process of computing all EFMs is based on 
the double description method [12 . The basic princi- 
pal of the double description method is to determine an 
initial set of solution modes which are then iteratively 
combined and added to the set of existing modes un- 



Rl 


-0.5 


0.5 


-0.5 0.5 


1 


0.5 


R2 


-1 


-1 


1 


1 





1 


R3 


1 

















R4 











1 


0.5 


0.5 


R5 

















1 


R6 














1 





R7f 





1 














R7b 








1 











R8 











1 








R9 














0.5 


0.5 


RIO 














1 





Rll 

















1 



TABLE II. Kernel matrix K of the extended stoichiometric 
matrix shown in Table [T] of the example network. 



til the complete set of modes is obtained. The solution 
modes are stored in the mode matrix R that contains one 
column for each elementary mode. Typically, the initial- 
ization of the mode matrix R is obtained by calculating 
the kernel K of the extended stoichiometric matrix S ex t- 
The kernel K of the extended stoichiometric matrix S ex t 
is defined by S ex tK = and is shown in Table |Tl| 

Next, the initial conversion to the binary representa- 
tion of the mode matrix R is performed. The final set 
of elementary modes of the extended network must only 
contain non-negative flux values, as the extended net- 
work contains only irreversible reactions. As pointed out 
by Gagneur and Klamt [1] using only irreversible reac- 
tions is of major importance, as all non-zero elements of a 
mode will be retained if a new mode is created by combin- 
ing this mode with other modes that have already been 
determined during the calculation procedure. Hence, all 
rows that contain only non-negative values can directly 
be transformed to the binary representation. For the sake 
of clarity we use the character / for binary 1 indicating 
a flux carrying reaction and the character n for binary 
indicating that no flux occurs. Usually, the initial solu- 
tion matrix R is sorted in a way that all rows containing 
only positive values are at the top. Table III shows the 



properly sorted mode matrix R containing numerical and 
binary values before the iteration process is started. 

Next, the iteration procedure is performed. Step by 
step each row that is still in numerical form is trans- 
formed to its binary representation. As shown in Ta- 
ble |III| the next reaction to be processed is R2. The 
double description method requires that all modes con- 
taining non-negative values at R2 are retained, whereas 
the modes with negative values are removed. Further- 
more, the method requires that all modes with negative 
values at R2 are combined with adjacent modes exhibit- 
ing a positive value at R2. Hence, the modes Ml and 
M2 are combined with M3, M4, and M6 resulting in six 
potential new modes. Two modes are adjacent if the 
binary part of the new mode is not a superset of any 
already existing modes - except the two parent modes. 
For the binary part, the combination of two adjacent 
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TABLE III. Initial mode matrix R. The first ten rows are al- 
ready transformed to the binary representation where /stands 
for binary 1 and indicates that the reaction carries a flux and 
n stands for binary and indicates that no flux is carried. 
Note that the order of the reactions has changed compared to 
Table U in order to maximize the number of rows that can be 
converted to binary form during the pre-iteration phase. 
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TABLE IV. Mode matrix R after the first iteration step con- 
verting reaction R2 from numerical to binary form. 



modes is a simple and fast bitwise OR operation of the in- 
volved modes. Combining the numerical part is achieved 
by a weighted subtraction of the two numerical vectors. 
The new numerical value v of row r is calculated by 

VneWr — {^posi^negr ^neg\^pos r ) I {Pposi ^negi ) i where 

v pos and v neg are the values of the positive and of the 
negative column at row r, respectively. The row index 
r runs from 1 to n, where row r = 1 is the row to be 
converted at current iteration step and n is the number 
of rows left to be converted. Applying these instructions 



to the initial mode matrix R given in Table III results in 
the new mode matrix shown in Table ITVl 

Applying the mode combination procedure again for 
the last row to be converted (Rl) results in the final mode 
matrix R as shown in Table |Vj Now, the matrix R con- 
tains only binary elements. Note that the performance of 
the described iteration procedure for 'real-life' networks 
can tremendously be increased if tree structures are uti- 
lized to store the binary representation of the modes [5] . 

Next, the futile 2-cycle mode (M6) that was created by 
splitting the reversible reaction R7r is removed. Then the 



TABLE V. Mode matrix R after the final iteration step con- 
taining only binary values. Mode M6 is the futile 2-cycle 
mode and can be removed. 





Rl 


R2 


R3 


R4 


R5 


R6 


R7r R8 


R9 


R10 


Rll 


EFM01 





1 





1 








1 


1 











EFM02 





1 





1 


1 





1 





1 





1 


EFM03 





1 





1 





1 


1 





1 


1 





EFM04 








1 


1 





1 


1 








1 





EFM05 


1 


1 





1 











1 











EFM06 


1 


1 





1 


1 











1 





1 


EFM07 


1 








1 





1 








1 


1 





EFM08 


1 








1 


1 





1 





1 





1 


EFM09 


1 








1 








1 


1 











EFM10 








1 


1 


1 











1 





1 


EFM11 








1 


1 











1 












TABLE VI. Binary representation of all elementary flux 
modes of the example network shown in Figure [I] 1 means 
that the reaction carries a flux and means the reaction car- 
ries no flux. Note that the futile two-cycle of the reversible 
reaction R7r has already been removed and the forward and 
backward irreversible reactions (i?7/and Rfb) have been com- 
bined to the reversible reaction R7rby a bitwise OR opera- 
tion. 



irreversible forward and backward reactions i?7/and Rib 
are combined by a simple bitwise OR operation in order 
to obtain the reversible reaction R7r again. The final set 

Recon- 



of modes in binary form is shown in Table VI 
verting the binary form to the numerical representation 
is achieved by using the fact that the reduced null space 
matrix N re d multiplied by the sought numerical mode has 
dimension 1 and is equal to zero. N rec i is a sub-matrix of 
the kernel K that only contains columns/reactions where 
the binary mode to be transformed carries a flux. Hence, 
only a homogeneous linear system has to be solved to ob- 
tain the 1-dimcnsional solution vector that represents the 
numerical form of a mode. The result of the reconversion 
of the binary modes for the simple example network is 
listed in Table Eni 

The binary approach combines several essential advan- 
tages: a) modes are stored in binary format which dra- 
matically reduces the memory usage, b) new modes are 
calculated from existing adjacent modes by using simple 
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TABLE VII. Numerical representation of all EFM of the ex- 
ample network shown in Figure [l] 



Of 




bitwise boolean functions which are very fast compared 
to numeric operations, and c) the bitwise boolean oper- 
ations used are 'exact', hence, numerical accuracy prob- 
lems are minimized. 



Gene regulation information 

TRN control the process of gene expression in cells. 
They determine how genes activate or repress certain 
fluxes. Hence, the gene regulatory information of a net- 
work imposes additional constraints on the reactions, 
and, as a consequence, reduces the solution space re- 
sulting in a lower number of biologically feasible EFMs. 
Typically, the gene regulation information is provided in 
form of boolean functions [llj . such as the NOT, OR, 
and AND operations. 

As illustrated in Figure [2] we assume that the gene reg- 
ulatory network of the example network shown in Figure 
[ljonly consists of a gene GR that activates reaction R7r 
and deactivates reactions R9. The function of gene GR 
can be transformed to a single boolean expression: R7r — 
NOT(i?P). This constraint means that the reaction R7r 
must not carry a flux when reaction R9 carries a flux and 
vice versa. A simple approach to get the reduced solution 
space is the application of this gene regulatory rule after 
all mathematically possible modes were calculated. Nat- 
urally, this method does not result in any performance 
improvement. However, if we consider the basic princi- 
ple of the binary approach described above, a significant 
speed up of the computation process can be obtained. 
The boolean operation R7r = NOT(i?9) implies that the 
rule is not obeyed if: a) R9 = 1 = / and R7r = 1 = /or 
b) R9 = = n and R7r = = n. The first expression 
is of particular interest, as it can be used to eliminate all 
modes from the solution matrix R - at any time of the 
iteration process - if R9 and R7r do carry a flux. This 
statement is true, as a) the considered mode itself dis- 
obeys the rule and b) all children modes generated from 
the considered mode by combination with other modes 



FIG. 2. Example network including the gene regulator net- 
work: R7r = NOT(i?3). 

will also disobey the rule. The latter property is owed 
to the fact that a flux value at a certain reaction will be 
retained by the bitwise OR operation for the rest of the 
computation procedure (see previous subsection for fur- 
ther details). Removing a mode as soon as possible is of 
essential importance, as this mode is hindered to create 
offspring modes that would have to be eliminated at a 
later stage. The second expression (if R9 = = n and 
R7r = = n) is of no use during the iteration process, 
as a no flux value of R9 or R7r can become a flux car- 
rying reaction in a child mode that is created in a later 
iteration step. Hence, removing a currently disobeying 
mode with R9 = and R7r = would result in the loss 
of children modes that obey the rule R7r = NOT(R9). 
However, the rule R7r = NOT{R9) for R9 = and R7r 
— can still be used to remove infeasible modes after 
finishing the computation of all binary modes 

The above considerations make clear that there are 
two types of rules: a) rules that can be applied during 
the iteration process and b) rules that can be applied 
during the post-processing step after finishing the mode 
calculation. 

Determining if a boolean rule Ro — B(Rl,...,Rn) qual- 
ifies for the iteration phase is simple. If the output re- 
action Ro of the rule is (does not carry a flux) when 
all input reactions Rl,...,Rn are 1 (carry a flux) then the 
rule can be used during the iteration phase. 

Special care must be taken for reversible reactions, as 
they are split and, hence, occur twice in the extended set 
of reactions. If either the forward or the backward reac- 
tion carries a flux then the original reaction is supposed 
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TABLE VIII. Mode matrix R after the first iteration step. 
The red font color highlights reactions involved in rule R7r 
= NOT(i?9) that carry a flux. Mode M5 (highlighted in grey 
background color) disobeys the rule and is removed from the 
matrix. 





Ml 


M2 


M3 M4 M5 


M6 


M7 


M8 




M9 




M10 


Mil 


R3 


n 


n 


n 


n 


n 


f 


f 


f 




n 




n 


n 


R4 


f 


f 


f 


f 


n 


f 


f 


f 




f 




f 


f 


R5 


n 


n 


f 


n 


n 


f 


n 


n 




f 




n 


n 


R6 


n 


f 


n 


n 


n 


n 


n 


f 




n 




f 


n 


R7f 


n 


n 


n 


f 


f 


n 


n 


n 




n 




n 


n 


R7b 


n 


n 


n 


n 


f 


n 


n 


f 




f 




f 


f 


R8 


f 


n 


n 


f 


n 


n 


f 


n 




n 




n 


f 


R9 


n 


1 


f 


n 


n 


f 


n 


f 




f 




f 


n 


R10 


n 


1 


n 


n 


n 


n 


n 


f 




n 




f 


n 


Rll 


n 


n 


1 


n 


n 


f 


n 


n 




f 




n 


n 


R2 


f 


n 


f 


n 


n 


n 


n 


n 




f 




f 


f 


Rl 


f 


1 


f 


f 


n 


n 


n 


n 




n 




n 


n 



TABLE IX. Mode matrix R after the final iteration step. 
Mode M8, M9, and M10 do not obey the iteration phase rule. 
Additionally, mode Ml and M7 disobey the post-processing 
rule. M5 is also removed, as it is the futile 2-cycle mode 
created by splitting the reversible reaction R7r into two irre- 
versible reaction. 



to be flux carrying when checked against a boolean rule. 

Applying these concepts to the example network with 
the gene regulatory rule R7r = NOT(R9) results in a 
mode matrix R after the first iteration step as shown 
in Table \VTU\ Table |V1H| highlights in red font color all 
reactions involved in rule R7r — NOT(i?<?) that carry a 
flux. It can be seen that mode M5 disobeys the rule and 
is removed from the matrix R. 

In the next step mode M5 does not exist and, hence, 
fewer adjacency tests have to be performed. Table IX 
shows the mode matrix R after the final iteration step. 
It can be seen that mode M8, M9, and M10 do not obey 
the iteration phase rule, as R9 and R7f or R7b carry a 
flux (highlighted by the red font color). Hence, M8, M9, 
and M10 can be removed during the iteration phase. 

Furthermore, Table IIXI illustrates that mode Ml and 
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0.5 1.0 0.0 0.5 1.0 0.0 0.0 0.0 0.5 0.0 1.0 

0.0 1.0 0.0 0.5 0.0 0.0 -0.5 0.5 0.0 0.0 0.0 

1.0 0.0 0.0 0.5 0.0 1.0 0.0 0.0 0.5 1.0 0.0 

1.0 0.0 0.0 1.0 0.0 0.0 1.0 1.0 0.0 0.0 0.0 

0.0 0.0 1.0 0.5 1.0 0.0 0.0 0.0 0.5 0.0 1.0 



TABLE X. Numerical representation of all EFMs of the ex- 
ample network shown in Figure [TJif the gene regulatory rule 
R7r = NOT(7?9) is applied. 



M7 disobey the post-processing rule, as R9, R7f, and R7b 
do not carry a flux value. Consequently, after removing 
the futile 2-cycle mode M5 the final mode matrix R only 
contains the five modes M2, M3, M4, M6, and Mil. Be- 
fore transforming the binary modes back to their numeri- 
cal form the split irreversible reactions i?7/and R7b must 
be combined to the reversible reaction R7r. The final set 
of feasible EFMs is listed in Table H 



Implementation 

We implemented our approach as an extension to the 
open source software efmtool. The mode elimination al- 
gorithm was realized by adding three Java packages to 
the original version of efmtool. The three packages con- 
tained ten new Java classes. These new classes are re- 
sponsible for handling the genetic rules and checking the 
modes against them. Two already existing Java classes 
were slightly enhanced in order to invoke the mode check. 
The boolean rules are provided by an additional input file 
using the command line argument -generule. The ex- 
tended version of efmtool was compiled by JDK 1.7.0. 
The implementation of the extension was performed as 
non-invasive as possible, which means that the perfor- 
mance gain might be even better if the new method is 
integrated to efmtool more thoroughly. The mode checks 
for the iteration phase were implemented using binary bit 
patterns where the patterns are created simply by setting 
the involved reactions (all input reactions and the out- 
put reaction) to 1 [13]. If a tested mode has every bits 
set that occurs in the binary bit pattern of a rule then 
the mode is eliminated. The mode check for the post- 
processing step was realized by utilizing a reverse polish 
notation approach that allows a simple and fast execu- 
tion of boolean functions with any values for the input 
reactions. The general sequence of operation of our ex- 
tended version of the binary double description method 
is shown in the supplementary data section. 



RESULTS AND DISCUSSION 

We tested our approach on the E. coli core model pro- 
vided by [TTJ [T3] . The model consists of 94 metabolites 
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Booleanly combined 


Required value of 


input reactions 


effected output reaction 


R_EX_glc_e = 1 


-> R_EX2_ac_e = 


R_EX_glc_e = 1 


-> R_ICL = 


R_EX_glu_L_e = 1 - 


-> R_GLUDy = 


R_EX_glu_L_e = 1 - 


-» R_GLUSy = 



TABLE XL The four boolean rules used by the introduced 
elimination algorithm to exclude biologically infeasible EFMs 
during the iteration phase. The 54 other rules can only be 
applied during the post-processing phase and are shown in 
the supplementary data section. 



and 95 reactions. 59 reactions are reversible. Gene reg- 
ulatory information for this model is provided by [H] in 
form of a gene-enzyme-reaction mapping. The mapping 
was checked for consistency and contradicting rules were 
removed from the provided data set. The final mapping 
used by our algorithm to exclude infeasible EFMs con- 
tained 58 boolean functions and is shown in Table 2 of the 
supplementary information section. Only four of these 58 
rules qualify for being used during the iteration phase of 
the EFM computation procedure and are listed in Table 

M 



Table XII compares a regular run without regulatory 
information and a run using the available gene regula- 
tion rules. Both runs were performed on a Linux Ubuntu 
11.04 computer with 2 Intel Xeon CPUs (6 cores each) 
and a total of 192 GB of RAM using 10 parallel threads. 
The table shows that after the iteration phase 226 mil- 
lion modes were obtained without regulatory informa- 
tion. Despite the fact that just four of the 58 boolean 
rules qualify for the iteration phase, applying the novel 
mode elimination approach resulted in only 76.7 million 
modes. This mode removal during the iteration phase 
caused a reduction of the iteration runtime from 30.9 to 
5.3 hours and a decrease of the maximum number of ad- 
jacent candidates from 2.2-10 15 to 2.6T0 14 . The elimina- 
tion of modes in the post-processing phase was even more 
pronounced, as only 2.1 million of the 76.7 million modes 
adhered to the boolean gene-enzyme-reaction mapping 
rules. Interestingly, even the post-processing runtime 
is reduced by our approach (from 3.2 to 1.8 hours) al- 
though 76.7 million binary modes were applied to the 
post-processing boolean rules. This behaviour is owed to 
the fact that only 2.6 million modes survive this post- 
processing filtering procedure which means that only a 
fraction of the modes were converted from the binary to 
the numerical representation and were written to disk. 
This huge reduction of the total number of modes (226.3 
million to 2.1 million) had a major influence on the re- 
quired harddisk space which was decreased from 251 GB 
to 2.3 GB if an uncompressed double precision text for- 
mat was used. Table [XXT] clearly shows that considering 
gene regulatory information in the computation process 
has a huge impact on the computational key properties 





w/o gene regulation 


with gene regulation 


No. of modes (iteration) 


226.3 ■ 10" 


76.7 ■ 10° 


No. of modes (post-processing) 


226.3 • 10° 


2.1 • 10 6 


Max. adjacent candidates 


2.2 • 10 15 


2.6 • 10 14 


Max. RAM usage 


153 GB 


73 GB 


Runtime (iteration) 


30.9 h 


5.3 h 


Runtime (post-processing) 


3.2 h 


1.8 h 


Runtime (total) 


34.1 h 


7.1 h 


Disk space 


251 GB 


2.3 GB 



TABLE XII. Comparison of EFM calculation with and with- 
out taking into account gene regulatory information. The re- 
quired disk space is given for a result file containing all modes 
in text format using double precision. The iteration runtime 
is the time spent creating the binary modes without pre- and 
post-processing and the line 'max. adjacent candidates' shows 
the maximum number of potentially occurring adjacent pairs. 



of the calculation of EFMs. 

In order to verify that the the presented extension of 
the efmtool computes the correct EFMs, an extra soft- 
ware tool was developed that applies the boolean rules 
to the complete and unfiltered set of EFMs. The two 
tools computed identical sets of EFMs ensuring that the 
efmtool extension produces the correct result. 

Tabic |XIII shows the development of the number of 
obtained modes as a function of finished iterations. In 
total, 21 iteration steps were performed in order to com- 
pute the complete set of EFMs. Up to iteration 9 not 
a single EFM was eliminated and the inclusion of gene 
regulatory information had no effect. The first removal 
occurred at iteration 10, where 3 modes were deleted. 
Although in total only 1.6 million modes were removed 
during the iteration phase, the final number of modes was 
reduced by 149.6 million modes. This huge reduction is 
a result of lost parent modes which otherwise could have 
spawned a multitude of new children. 

In order to find an initial value of the mode matrix R 
the kernel of the extended stoichiometric matrix is com- 
puted. Before the iteration phase is started the reactions 
are sorted. This is done to put all reactions with only 
positive values to the top which results in the maximum 
number of reactions that can be transformed to the bi- 
nary form before the iteration procedure is even started. 
Several approaches can be applied to sort the reactions 
that also contain negative values, e.g. taking no spe- 
cial measures (random order) or ordering by increasing 
potential combinations (number of negative values times 
number of positive values). As the iteration phase rules 
can only be applied if the involved reactions are already 
converted to the binary representation, it seems benefi- 
cial to convert all reactions that are involved in iteration 
phase rules to the binary form as soon as possible. This 
concept requires that all these reactions are moved to the 
top of the set of rows containing negative values. Note 
that this approach was not implemented in the developed 
extension but could result in an additional performance 
gain if realized. 
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Iteration 


No. of removed 


No. of modes incl. 


No. of modes w/o 


No. 


infeasible modes 


gene regulation 


gene regulation 


9 





97 


97 


10 


3 


205 


208 


11 





285 


288 


12 





454 


457 


13 


o 


456 


459 


14 


o 


751 


755 


15 


o 


849 


952 


16 


604 


2,113 


3,223 


17 


3 


6,463 


9,454 


18 


850 


17,154 


28,114 


19 


2,168 


27,468 


48,388 


20 





27,468 


48,388 


21 


1,597 


57,180 


112,180 


22 


1,717 


244,858 


486,847 


23 


81 


224,537 


444,371 


24 


93 933 


519 853 


1.243,347 


25 


109,042 


1 701 029 


4 566 570 


26 


295,410 


2,832,654 


8,012,612 


27 


31,045 


4,505,295 


12,790,524 


28 


250,704 


12,895,654 


37,465,244 


29 


247,738 


26,365,168 


77,934,795 


30 


59,107 


32,421,087 


94,929,161 


31 


591,718 


76,690,502 


226,269,046 


sum 


1,685,720 







TABLE XIII. Comparison of the number of EFMs as a func- 
tion of the iteration step for simulations with gene regulatory 
information and without gene regulatory information. 



CONCLUSION 

We implemented a novel approach to speed up the 
computation of the complete set of EFMs of a metabolic 
network by extending the open source program efmtool 
written by Marco Terzer. Our extension allows the con- 
sideration of gene-enzyme-reaction mappings in the pro- 
cess of EFM computation. Biologically infeasible flux 
modes are constantly eliminated during the calculation 
process. By implementing an early stage exclusion of 
modes a considerable reduction of computational costs 
was achieved which pushes the maximum size of calcula- 
ble networks to new limits. We think that our approach 
is another step to the final goal of studying genome-scale 
metabolic networks by elementary flux modes. 
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